Modeling and simulation of synthetic nanopores using MsSimPore 
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Abstract 

MsSimPore is a novel finite element software package for the simulation of ion flow through 
nanopores, which are radially symmetric or whose lateral extension is negligible to the longi- 
tudinal one. Its graphical user interface allows for the variation of experimental parameters 
and the subsequent visualization of the calculated results. MsSimPore offers different solver 
options, either simulations of the concentration profiles for single applied voltages or the 
calculation of current-voltage and rectification curves. 

In this paper we discuss the basic features of MsSimPore and the underlying mathematical 
modeling. Furthermore we compare the simulations with current voltage curves, obtained 
from conical and cigar shaped nanopores, under different experimental conditions. 

1 Introduction 

Nanopores are nanoscale channels in synthetic materials such as silicon nitride, graphene or 
polymers, e.g. polyethylene terephthalate [U [2], [3J SJ, [5j (6J El [8] . They can be made in a variety 
of lengths, diameters and shapes, only limited by the thickness and robustness of the membrane 
material. The transport properties of the pores can be tuned and modified by the surface 
charge of the material and by chemical modifications of the channel surface, e.g. via attaching 
large biomolecules. Due to their versatile and robust behavior, nanopores have emerged as 
promising tools for regulating the transport of charged particles, sensing single molecules or 
DNA sequencing in the last years P fTUl fTTl IT21 [T3] . 

The Poisson-Nernst-Planck (PNP) equations have been used successfully for the simulation of 
ion flux through solid state nanopores |14|, Ha] [TE[ I17j . Cervera and co-workers presented first 
ID simulations in good agreement with experimental data for conically shaped nanopores |14j . 
Furthermore they discussed transport and rectification properties for cigar shaped nanopores. 
Here the results were solely based on simulations in ID of a reduced PNP system, where the 
Poisson equation for the potential was replaced by the charge neutrality condition [15J. In 
both cases simulations were performed in the pore region only, hence they assumed that the 
influence of the access resistance can be neglected |18] . This assumption results in the Donnan 



equilibrium values for the ion concentrations and the local electroneutrality to determine the 
values of the potential. More recently different generalizations of the PNP equations, which 
include the ionic size, have been proposed in the literature. This can be accomplished via 
an additional potential in the transport equation[19j. A different approach includes the ion 
size already in the derivation, which leads to nonlinear mobilities and diffusivities |20 | 116 1 f2Tj . 
Another generalization of the PNP equations in 2D, which includes chemical reactions to describe 
the formation of nanoprecipitates in nanopores, has been proposed by Wolfram and co-workers 
[22] . 

The software package MsPoreSim is based on a one-dimensional reduction of the classical PNP 
equations. It allows accurate simulations of the full PNP system of nanopores, which are either 
radially symmetric or whose transversal dimension can be neglected, with several ionic species 
present in the bath. The computational domain includes the bath and pore region, giving a more 
accurate description of the experimental setup. MsSimPore is based on a robust finite element 
solver for the stationary PNP equations, which allows simulations of versatile pores of various 
geometries and external conditions. Its graphical user interface is easy to use and automatically 
generates graphical images of the concentration and voltage profiles, as well as current-voltage 
(IV) and rectification curves. The fast and stable simulations allow the use of MsSimPore for 
parameter fitting problems, e.g. to fit IV curves by changing different parameters, like the 
surface charge or the shape of the pore. 

In this paper we present the underlying mathematical model and numerical algorithms used. 
Furthermore we compare the numerical simulations with various IV curves for differently shaped 
nanopores under various experimental conditions. In particular we are interested in the effects 
of the bath concentration, pH value and temperature on the surface charge of the pore wall. 
Shielding of the surface charges on the pore walls, due to high bath concentrations, has been 
observed in different experimental settings |23t 124] . Also the surface charge decreases due to 
protonation at lower pH values [25]. This behavior is confirmed when fitting the experimental 
data, i.e. to match the given IV curves, the surface charge for higher bath concentrations and 
lower pH values needs to be decreased. 



2 Methods 
2.1 Experimental 

The nanopores were fabricated by irradiating 12 thick polyethylene terephthalate (PET) foils 
with exactly one single heavy ion of GeV energy. Each projectile produces a so-called ion track 
consisting of damaged material of few nm in diameter [26]. The track in the foil is converted 
into an open channel by chemical etching [271 128j . For this, the foil is mounted between two 
chambers of a custom-made conductivity cell, with one chamber being filled with 9 M NaOH and 
the other one filled with a neutralizing solution [29] . Given by the high NaOH concentration, 
dissolving the track from one side, conical nanochannels are created. Replacing the etchant by 
an electrolyte, current- volt age curves were recorded in the same cell using Ag/AgCl electrodes 
(chloridated Ag wires), a Keithley 6487 picoammeter /voltage source, and various KC1 solutions 
(stock solution of 1 M KC1, lower concentrations prepared by dilution). The electrode at the 
small opening of the pore was grounded, while the other electrode, placed in the cell chamber 
with the large opening of the pore, was used to apply a given transmembrane potential with 
respect to the ground electrode. 
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2.2 Modeling and simulation 

MsSimPore has been developed for conical and cigar shaped nanopores, as described in the 
previous section. The shape of the pores exhibit radial symmetry, therefore we can use a ID 
reduction of the PNP equations [30J. The main driving forces in the classical PNP equations 
are diffusion and electrostatic interactions with other ions, as well as the surface charges on the 
pore walls. 

In analogy with the experimental setup, MsSimPore simulates a single pore of length L separat- 
ing two electrolyte solutions. The electrolyte may have different concentration on the left and 
right-hand-side. The small and large opening radius of the conical pore is denoted by r s and 
ri, respectively. The computational domain considers SI = [0, 5L], i.e. the pore separates two 
electrolyte solutions, each container has a size of 2L (this is usually sufficient for equilibration) 
attached to the right and left. The electric potential is given by V = V(x), where x corresponds 
to the position in space. The small opening of the pore is positioned at x = 2L, the large 
at x = 3L. The concentration of each ionic species present in the electrolyte is pi = Pi(x), 
i = 1, . . . m. The reduced PNP model reads as: 

- dxv{eA{x)VV) = eA(x) ^ Ztpi + dA(x) cr(x) (la) 

i 

= div(A(x)D l (x)(V Pi + Zi-^-piVV)), (lb) 

k B T 

where A = A{x) describes the cross-section of the pore and dA = dA(x) its circumference. 
Here e denotes the dielectric coefficient, e the elementary charge, fee the Boltzmann constant, 
T the temperature, z% and Di the valence and diffusivity of each ionic species, respectively. The 
function a = a(x) corresponds to the surface charge inside the pore. 
The area function is given by 

oi - a 2 exp(-(L//i)") - (r s - n ) exp(-(x/L)"(L/fr)") 
(X) l-exp(-(L/h) n ) 

inside the pore, and takes a large value in the bath regions [16]. The constants a\ and a 2 are 
chosen such that A(x) = r^ir at x = 2L and A(x) = rfn at x = 3L. The ratio (L/h) and the 
parameter n determine the curved shape of the pore. If {L/h) — > 0, the area function A = A(x) 
corresponds to the linear interpolation between circles of radius r s and ri, modeling a conical 
pore. For large ratios the shape of the pore is more curved, looking like a cigar, see also [T] 

The bath concentrations of each ionic species are modeled by Dirichlet boundary conditions, 
hence we have pi(x) = pi at x = and x = 5L for each i = 1, . . . m, where pi denotes the molar 
concentration of each ionic species in the bath. Also the applied voltage V app i is modeled via a 
Dirichlet boundary condition, i.e. V = V app i at x = and x = 5L. 

Two well known reformulations of [T] can be found in the literature, either based on the Slotboom 
variables or the entropy variables (also known as quasi Fermi potentials in the semiconductor 
community). The Slotboom variables, Uj = piexp(cziV) guarantee positive concentrations, but 
the exponentials can cause overflow problems for large applied voltages V app i . This problem can 
be avoided using entropy variables cpi = log pi + cziV. Then [T] reads as 

-A 2 div(eA(x)W) = kA(x) ^ z iPl + dA{x)cr{x) (2a) 

i 

= div(Di{x)A(x)(piV(pi)). (2b) 

Here A 2 = j^r^ , k = and c = eV/ksT are scaling parameters, and V,A,... denote typi- 
cal values of the physical constants, see [TJ This non-dimensionalization allows stable and unit 
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Figure 1: Sketch of nanochannel separating two electrolyte containers. Possible area functions 
are shown (not in scale), where the pore shapes correspond to different values of the parameter 
(L/h). The curved shape of the small pore opening depends on the ratio of (L/h), the smaller 
the ratio the more linear the interpolation between the small and the large opening radius. 



Meaning 


Value 


Unit 


Boltzmann constant fee 


1.3806504 x 10~ 23 


J/K 


Vacuum permittivity eo 


8.854187817 x 10~ 12 


F/m 


Relative permittivity e r 


78.4 




Elementary charge e 


1.602176 x 10~ 19 


C 


Temperature T 


293.16 


K 


Typical length L 


1 


nm 


Typical concentration c 


3.7037 x 10 25 


N/l 


Typical voltage V 


100 


mV 



Table 1: Parameters for computation 



independent simulations. 

We solve system [2] on = [0, 5L], where [0,2L] and [3L,5L] correspond to the left and right 
bath respectively, [2L, 3L] is the pore region. [2] is discretized using a hybrid discontinuous 
Galerkin method with upwind stabilization[31j. This stabilization ensures stability of the nu- 
merical scheme for large applied voltages. The discrete nonlinear problem is solved by Newton's 
method. The calculation of the Newton update is based on the non-symmetric solver MUMPS 
[321 [33]. 

MsSimPore has been implemented within the finite element framework of Netgen/NgSolve [34J 
and uses Qt|35| for its graphical interface. It allows the simulation of conical and cigar shaped 
nanopores for up to six species present in the bath. The graphical user interface distinguishes 
between pore related input parameters (e.g. small and large opening radius r s and r;, surface 
charge <r,...) and the experimental conditions (e.g. number of species present in the bath and 
their respective concentrations pi, applied voltage V app i, temperature....). We assume that the 
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Figure 2: Simulated (points) and experimental IV curves (solid lines) for a conical nanopore 
with radii r s = 9 nm and r\ = 450 nm for symmetric KC1 bath concentrations 0.01, 0.1, 0.2, 0.5 
and 1 M 



surface charge a is constant inside the pore, which seems a reasonable assumption for the pores 
considered. The ion specific parameters, like the valence and the diffusion coefficient, are stored 
in a small data base, which can be modified and extended by the user. MsSimPore offers sev- 
eral solver options such as simulations for one particularly chosen applied voltage as well as 
the calculation of current-voltage (IV) and rectification curves. The solver output variables are 
displayed in the graphical user interface and stored in a neutral format for subsequent personal 
use. 

MsSimPore uses an adaptive mesh, i.e. a coarse discretization of the computational domain O in 
the bath regions (using a mesh size of h = 50 nm), which we refine (as small as h m i n = 0.1 nm) 
around the narrow tip to resolve the fine features correctly. This automatic refinement reduces 
the computational costs and allows faster simulations. MsSimPore is available for download at 
the University of Minister [36J. 

3 Examples and discussion 

We illustrate and compare the results of modeling ion current with MsSimPore, with exper- 
imental data recorded for nanopores of different geometry and under different experimental 
conditions. Here we assume that all pores were radially symmetric and use IV curves recorded 
under symmetric bath concentrations of KC1. The diffusion coefficients for potassium and chlo- 
ride are set to Dk = 1-95 x 10 -9 m 2 /s and Dei = 2.03 x 10~ 9 m 2 /s. 

We consider a conically shaped nanopore with opening radii r s = 9 nm and r\ = 450 nm, 
which separates two solutions of KC1. [2] shows the IV values (symbols) generated by MsSimPore 
and the corresponding experimental measurements (solid line) for different bath concentrations. 
Best agreement with the experimental data is achieved, if we choose a = — 0.4e/nm 2 for 0.01 
M bath concentrations and a = — 0.05e/nm 2 for 1 M bath concentrations, instead of the ini- 
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Applied potential [V] 

Figure 3: Simulated (points) and experimental IV curves (solid line) for a conical nanopore 
with radii r s = 6 nm and 77 = 812 nm. The first three fits (blue, red and green) correspond 
to symmetric KC1 concentrations 1, 0.1 and 0.01 at pH 8, the last two (grey and orange) to 
symmetric KC1 concentrations of 0.1 and 0.01M KC1 at pH 5.5 



tially assumed value of a = — 1 e/nm 2 . Different experiments confirm that nanopores loose their 
rectification property for large bath concentrations, showing an almost linear behavior for large 
molar concentrations. This can be explained by the screening of the surface charges occurring 
over a small distance at high salt concentrations, thus reducing the effect of surface charges on 
the ionic concentrations in the pore. 

In [3] we consider a conical pore with radii r s = 6 nm and ri = 812 nm, where IV curves have 
been measured at different pH values. Again we obtain better fits if we decrease the surface 
charge for larger concentrations, i.e. a = —0.5 e/nm 2 for 0.01 M KC1 and a = —0.1 e/nm 2 for 
1 M KC1, see[3j At a lower pH value these values decrease even further to a = —0.1 e/nm 2 for 
0.01 M KC1 and a = —0.05 e/nm 2 for 0.1 M KC1. This decrease is caused by the protonation 
of the negatively charged carboxyl groups of the polymer wall. At even lower pH values the 
surface charge would further diminish. 

Finally we illustrate and compare the results of the solver for a cigar shaped pores, depicted in 
[4} Here only the large opening radius is known to be about ~ 180 nm. This value is based 
on etching multiple pore membranes, but at this pore size the uncertainty of the big opening 
diameter is about 10 %. All other parameters, like the small radius r s and the ratio determining 
the curvyness L/h, have been fitted using the experimental IV curves. As the number of un- 
known parameters increases, the fitting process plays a more dominant role. While for conical 
shaped nanopores, the fitting can be done by hand, the process is less clear for cigar shaped 
pores. Because of too many parameters being involved, it is not so clear, which parameter set 
fits the experimental data better. 

The experimental data shows strong rectification and compared (to conical pores) large currents 
for small bath concentrations. This indicates that the pore has probably a curved opening and 
that the r s is in the nanometer range [15j. [5] shows two different fits using similar parameter 
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Figure 4: Scanning electron microscope (SEM) image of cross-section of cigar shaped nanopore 
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Figure 5: Two fits for cigar shaped nanopore: simulated (points) and experimental IV curves 
(solid line). 



sets. In the first case the parameters are r s = 2nm, ri = 187nm and L/h = 100, while in the 
second fit the values are r s = lnm, r\ = 200nm and L/h = 100. The corresponding surface 
charges were a = —2, —0.2, and —0.05 e/nm 2 inland a = —2.5, —0.3, and —0.05 e/nm 2 in[3]for 
0.01,0.1, and 1 M KCI bath concentrations, respectively. In [3] the simulations are in excellent 
agreement with the experimental data for 0.01 and 0.1 M KCI, but deviate for 1 M KCI. If 
we decrease the large radius to 187 nm and increase the small radius to 2 nm the fit for 1 M 
KCI improves, but the curves for 0.1 and 0.01 start to deviate. This example demonstrates 
the complexity and difficulty of fitting data sets with several unknown parameters. To obtain 
more reliable results, advanced mathematical techniques for parameter identification problems, 
developed within the framework of inverse problems, have to be used. The latter experiment 
indicates an important future direction of possible mathematical research and software devel- 
opment, namely to determine unknown structural parameters in a reliable and automatic way. 



4 Conclusion 

We presented the software package MsSimPore developed for the simulation of conical and 
cigar shaped nanopores. The underlying mathematical model is a one-dimensional reduction 
of the steady state PNP equations. MsSimPore allows reliable, stable and fast calculations of 



concentration and voltage profiles as well as IV and rectification curves. We confirmed the 
computational results using IV curves of differently shaped pores under various experimental 
conditions, see [2] and [3j 

Further developments of MsSimPore will mainly focus on identification or design problems, where 
physical or structural parameters of nanopores are calculated from indirect measurements like 
IV or rectification curves. The presented examples, see [5j already indicated two identification 
problems, which arise in nanopore simulations, i.e. the identification of the surface charge 
a = <j{x) and the small opening radius r s . Both parameters are difficult to access experimentally. 
The surface charge of nanopores was only deduced indirectly and is estimated to be — 1 e /nm 2 . 
The shape of the pore can be analyzed by the replica technique [371 EE] , where the pore is filled 
galvanically with a metal. Subsequent dissolution of the polymer matrix reveals the replicated 
pore geometry. During this process, the narrow tip (which influences the behavior of the pore 
significantly) often breaks and the precise shape of the small pore opening can not be determined. 
Its size and geometrical shape could be identified from indirect measurements like IV curves 
using methods from inverse problems. Burger et al. presented first results on inverse problems 
in ion channels |39j and also a comprehensive overview and lookout on problems and methods 
for identification problems in synthetic and biological pores [4Q|. |4"T] . This direction of research 
is a promising application of advanced mathematical techniques for the automated, stable and 
reliable identification of structural parameters in nanopore experiments. 
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